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ABSTRACT 


This  report  presents  a  hybrid  method  for  simulating  sequences  of  correlated  Gamma 
random  variables  for  modelling  sea  clutter,  using  a  combination  of  linear  and/or  non¬ 
linear  transforms.  Depending  on  the  shape  parameter,  this  method  minimises  the 
use  of  non-linear  transformations.  Mathematically  the  method  is  simpler  than  its 
counterpart  methods  which  leads  to  a  quicker  simulation  run  time.  Two  memoryless 
non-linear  transform  (MNLT)  approaches  are  also  studied  with  comparative  results 
showing  that  the  hybrid  approach  is  more  computationally  efficient  and  slightly  more 
accurate  for  low  shape  parameters.  The  drawback  of  the  proposed  method  is,  however, 
that  it  can  only  handle  positive  correlations  whilst  the  two  MNLT  methods  are  capable 
of  handling  both  positive  and  negative  correlations. 
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Generating  Correlated  Gamma  Sequences  for  Sea-Clutter  Simulation 

Executive  Summary 

To  support  the  Generic  Phased  Array  Radar  Modelling  (GPARM)  simulation  program  of  the  De¬ 
fence  Science  and  Technology  Organisation  (DSTO),  which  in  turn  supports  many  Australian 
defence  programs  including  SEA1448  (ANZAC  ASMD),  AIR7000  (future  maritime  patrol  and 
response  capability)  and  AIR5077  (Wedgetail),  this  report  proposes  a  hybrid  method  and  exam¬ 
ines  other  existing  methods  for  simulating  sequences  of  correlated  random  variables  with  a  Gamma 
distribution. 

Radar  sea-clutter  is  a  dominant  undesired  signal  seriously  affecting  radar  performance  over 
the  sea  surface.  It  is  now  widely  accepted  that  in  most  scenarios,  radar  sea-clutter  can  be  modelled 
as  a  compound  non-Gaussian  random  process.  One  of  the  most  popular  is  the  compound  Re¬ 
distribution  which  consists  of  two  parts:  a  fast- varying  component  representing  sea-clutter  speckle 
which  is  modelled  as  a  complex  Gaussian  with  zero  mean  and  unit  variance  and  a  slowly-varying 
component  which  is  Gamma  distributed  and  represents  the  underlying  sea-clutter  intensity.  These 
two  components  are  assumed  to  be  mutually  independent. 

Depending  on  radar  parameters  and  sea  surface  conditions,  each  component  of  the  received 
sea-clutter  may  be  correlated  and  appropriate  correlation  models  should  be  included  in  the  simu¬ 
lation.  For  K-distributed  sea-clutter,  methods  for  simulating  both  correlated  Gaussian  and  Gamma 
processes  are  required.  Simulation  of  the  former  process  is  straightforward  and  can  be  realised 
by  a  linear  transform  using  either  spherically  invariant  random  processes  (SIRP)  or  Fourier  syn¬ 
thesis.  The  advantage  of  the  linear  transform  is  that  the  desired  correlation  properties  are  easily 
maintained.  However,  simulation  of  the  correlated  Gamma  distribution  is  more  difficult  and  may 
require  application  of  the  so-called  memoryless  non-linear  transform  (MNFT)  to  generate  the  de¬ 
sired  correlation. 

This  report  proposes  a  hybrid  method  for  simulating  sequences  of  correlated  Gamma  random 
variables.  The  approach  depends  on  both  the  desired  shape  parameter  and  in  some  cases  the 
correlation  coefficient.  For  most  distributions  with  a  shape  parameter  greater  than  0.5,  the  method 
only  requires  linear  transforms  to  generate  the  desired  Gamma  correlation.  However  in  other  cases, 
such  as  when  the  shape  parameter  is  less  than  0.5,  the  method  requires  the  MNFT  to  achieve  the 
desired  correlation.  The  MFNT  technique  proposed  in  this  report  is  different  from  its  counterparts. 
Compared  to  the  other  methods,  the  proposed  one  is  mathematically  simpler  making  its  numerical 
implementation  easier  and  more  computationally  efficient.  The  drawback  is,  however,  that  it  can 
only  handle  positive  correlation  coefficients. 

Two  implementations  of  the  MNFT  are  also  examined  and  evaluated  in  this  report.  The  first 
directly  implements  the  numerical  integration  described  by  Tough  and  Ward,  while  the  second  is 
a  polynomial  auto-correlation  method  extended  from  the  former  by  Weinberg  and  Gunn.  Results 
for  the  former  have  a  high  level  of  accuracy  and  speed  due  to  the  optimised  code.  The  latter  also 
performs  reasonably  well  except  for  small  shape  parameters  less  than  0.3.  Whilst  these  methods 
are  much  more  complex  then  the  hybrid  method,  they  are  able  to  handle  negative  correlation  coef¬ 
ficients.  A  comparison  between  methods  shows  that  the  hybrid  method  has  a  reduced  simulation 
run  time  and  provides  slightly  higher  accuracy  in  generating  the  desired  correlated  sequences, 
particularly  for  small  shape  values. 
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Another  contribution  in  this  report  is  the  provision  of  criteria  for  realisable  correlation  func¬ 
tions.  It  is  shown  that  not  every  correlation  function  can  be  simulated  but  only  those  that  have  a 
positive  semi-definite  or  positive  definite  covariance  matrix. 

This  report  serves  as  a  detailed  reference  for  computer  programmers  and  research  scientists 
who  want  to  implement  correlated  Gamma  processes  to  simulate  distributions  such  as  the  Re¬ 
distribution. 
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1  Introduction 


The  compound  K-distribution  is  often  used  to  model  radar  sea-clutter.  It  consists  of  a  fast- varying 
component  representing  sea-clutter  speckle  that  is  modelled  as  a  Gaussian  process  with  zero  mean 
and  unit  variance  and  a  slowly- varying  component  representing  fluctuation  of  the  underlying  sea- 
clutter  intensity  that  is  Gamma  distributed. 

Depending  on  radar  parameters  and  sea  surface  conditions,  the  received  sea-clutter  may  be  cor¬ 
related,  with  each  of  the  two  components  making  contributions  to  the  correlation.  Using  a  pulsed 
Doppler  radar,  the  received  sea-clutter  data  in  a  coherent  processing  interval  (CPI)  is  in  general  a 
two-dimensional  data- sef]  (time-range  or  Doppler-range).  The  temporal  correlation  refers  to  the 
correlation  between  pulses  separated  in  time,  while  the  spatial  correlation  refers  to  the  correlation 
between  range  bins.  The  definition  and  calculation  of  both  the  temporal  and  spatial  correlations  is 
given  in  Appendix  |A| 

Methods  for  simulating  correlated  K-distributed  sea-clutter  are  a  highly  relevant  topic  for  the 
radar  community.  Since  the  fast-varying  component  is  Gaussian  distributed,  simulating  its  cor¬ 
relation  is  relatively  easy  using  the  well-established  spherically  invariant  random  process  (SIRP) 
[Rangaswamy,  Weiner  &  Ozturk  1993,  Antipov  1998].  On  the  contrary,  the  simulation  of  cor¬ 
related  non-Gaussian  processes  is  not  straightforward.  Methods  for  simulating  these  stochastic 
processes  can  be  classified  into  two  categories.  The  first  uses  the  SIRP  formulation  [Rangaswamy, 
Weiner  &  Ozturk  1993]  and  requires  the  processes  to  be  linearly  transformable  from  the  Gaussian 
process.  To  implement,  a  linear  transform  relationship  is  first  established  between  the  desired  and 
Gaussian  distributions  and  then  the  correlation  of  the  desired  distribution  is  mapped  onto  the  cor¬ 
relation  of  the  Gaussian  process.  The  correlated  Gaussian  process  is  then  realised  through  SIRP 
before  being  transformed  back  to  the  desired  distribution.  Since  the  linear  transform  maintains  the 
correlation  properties,  SIRP  can  achieve  a  high  level  of  accuracy.  However,  only  a  few  random 
stochastic  processes  can  be  linearly  connected  to  the  Gaussian  process. 

The  generation  of  correlated  Gamma  random  fields  via  SIRP  theory  is  examined  in  [Conte 
et  al.  1991,  Armstrong  &  Griffiths  1991].  In  these  papers,  the  Gamma  random  process  is  ex¬ 
pressed  as  a  sum  of  Gaussian  processes  when  the  Gamma  shape  parameter  v  has  special  values. 
The  method  of  [Conte  et  al.  1991]  shows  that  the  Gamma  distribution,  in  the  range  of  0  <  v  <  2, 
closely  resembles  a  Gaussian  process  and  hence  the  correlated  Gamma  with  a  shape  parame¬ 
ter  in  that  range  can  be  produced  via  SIRP.  On  the  other  hand,  the  method  of  [Armstrong  & 
Griffiths  1991]  only  considers  a  special  correlation  case  where  the  correlation  follows  a  geometric 
progression  (i.e.,  the  correlation  coefficients  form  a  geometric  sequence)  for  v ,  having  values  of 
v  —  0.5n  where  n  is  a  natural  number. 

The  second  category  of  methods  for  simulating  correlated  random  stochastic  processes  is  the 
memoryless  non-linear  transformation  (MNLT).  This  method  is  examined  extensively  in  [Tough 
&  Ward  1999]  for  simulating  different  correlated  non-Gaussian  variates  including  Gamma.  The 
difficulty  of  the  MNLT  lies  in  the  non-linear  mapping  of  the  auto-correlation  of  the  output  (the 
desired  non-Gaussian  process)  to  the  auto-correlation  of  the  input  (the  Gaussian  process). 

This  report  proposes  a  hybrid  method  for  simulating  correlated  Gamma  random  variables.  The 
approach  depends  on  both  the  desired  shape  parameter  and  in  some  cases  the  correlation  coeffi¬ 
cient.  For  most  distributions  with  a  shape  parameter  greater  than  0.5,  the  method  only  requires 

'it  can  be  three-dimensional  if  radar  has  a  multi-channel  receiver. 
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linear  transforms  of  SIRP  to  generate  the  desired  Gamma  correlation.  However  in  other  cases, 
such  as  when  the  shape  parameter  is  less  than  0.5,  the  method  requires  the  MNLT  to  achieve  the 
desired  correlation. 

Two  implementations  of  the  MNLT  are  also  examined  and  evaluated  in  this  report.  The  first 
directly  implements  the  approximation  of  the  non-linear  mapping  described  in  [Tough  &  Ward 
1999].  This  method  is  referred  to  as  the  MNLT  with  numerical  integration.  The  second  method  is 
described  in  [Weinberg  &  Gunn  2011a]  and  extends  this  method  with  a  polynomial  approximation 
to  allow  near  optimal  control  between  the  input  and  output  processes.  This  method  is  referred  to 
as  the  MNLT  with  polynomial  approximation. 

The  report  is  organised  as  follows.  In  Section  [2}  the  proposed  hybrid  method  for  simulating 
correlated  Gamma  random  variables  is  described.  The  method  extends  those  of  [Conte  et  al.  1991, 
Armstrong  &  Griffiths  1991]  to  enable  an  arbitrary  shape  parameter  to  be  obtained  through  SIRP 
for  most  cases.  Details  of  the  two  MNLT  methods  are  then  briefly  described  in  Section  [3]  with 
evaluation  and  comparison  of  the  three  methods  presented  in  Section  |4j 
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2  The  Hybrid  Method 


A  random  variable  Z  is  said  to  have  a  Gamma  distribution  if  it  has  a  probability  distribution 
function  (PDF)  of, 


Pz(z)  = 

where  7(1/)  =  {y  —  1)!  is  the  Gamma  function.  The  above  Gamma  distribution  is  commonly 
denoted  by  Z  ^  r(z/,  9 ),  where  v  is  the  shape  parameter  and  9  the  scale  parameter.  It  is  known 
that  the  Gamma  distribution  has  the  following  property: 

Lemma  7:  If  Zk  ~  r(z/^,  9)  for  k  =  1, . . .  K  (i.e.,  all  have  the  same  scale  parameter,  but  may 
have  different  shape  parameters),  and  all  Z k  are  mutually  independent,  then  J2k= 1  ~  T  (z 7  9), 

where  v  =  Ylk= 1  uk- 

Lemma  1  only  requires  that  Z \  are  mutually  independent  and  not  whether  sequences  of  Z \ 
are  auto-correlated.  Consider  samples  of  the  random  process  Zk  taken  at  intervals  t  =  nAt 
where  1/A^  is  the  sampling  rate.  A  discrete  sequence  is  defined  by  Zk[n]  ~  T  {y,  9)  where  n 
is  a  natural  number.  The  following  two  properties  reveal  the  relationship  between  Gamma  and 
Gaussian  processes  and  thus  SIRP  can  be  established  from  Gaussian  to  Gamma. 

Property  1:  If  X  is  a  real  Gaussian  variable,  X  ~  N( 0, 0/2),  then  its  intensity  Z  =  \X\2  has  a 
Gamma  distribution  of  Z  ^  r(0.5,  9). 

Property  2:  If  X  is  a  complex  Gaussian  variable,  X  ~  CN( 0,  6 ),  then  its  intensity  Z  has  a 
Gamma  distribution  of  Z  r  (M). 


7  (v)6' 

0 


-,zv~1exp 


(-!)• 

2:  <  0, 


(1) 


Strictly  speaking,  sea-clutter  in  general  is  not  a  wide-sense  stationary  (WSS)  process.  However, 
the  correlation  of  the  Gamma  can  be  assumed  stationary.  The  goal  of  the  hybrid  method  is  to 
generate  a  sequence  z[n\  that  has  a  T  (y,  9)  distribution  and  a  correlation  coefficient  pk  given  by, 


Pk 


E{z[n\z[n  +  k ]}  —  z2 
var  (z) 


k  =  0, 1, . . .  and  p_k  =  Pk- 


(2) 


where  var (z)  is  the  variance  of  z.  The  relationship  between  the  covariance  E{(z[n\  —  z)(z[n  + 
k\  —  z)}  and  the  correlation  coefficient  pk  is  governed  by  (JT|).  Since  z[n\  is  Gamma  distributed, 
as  Z  rsj  T(y,9),  z  =  v9 ,  z2  =  v92  +  v292  and  var(^)  =  v9 2,  then  po  =  1  and  \pk\  <  1  for 
k  >  0.  Note  that  if  the  correlated  Gamma  is  formed  by  linear  correlated  Gaussian  processes,  the 
lowest  bound  of  pk  achievable  is  inherently  0  which  means  E{z[n]z[n  +  k}}  >  E2{z[n]}.  This 
will  become  clear  in  the  next  subsectioi0 


The  approach  of  the  hybrid  method  utilises  these  properties  to  efficiently  produce  the  corre¬ 
lated  random  sequence.  The  exact  approach  depends  on  both  the  desired  shape  parameter  and 
in  some  cases  the  correlation  coefficient.  For  most  distributions  with  a  shape  parameter  greater 
than  0.5,  the  method  only  requires  linear  transforms  of  SIRP  or  Fourier  synthesis  to  generate  the 
desired  correlation.  However  in  other  cases,  such  as  when  the  shape  parameter  is  less  than  0.5,  the 
method  requires  the  MNLT  to  achieve  the  desired  correlation.  The  following  three  cases  describe 
this  in  more  detail. 

2However,  if  correlated  Gamma  is  formed  by  the  MNLT  methods  discussed  in  Section [3]  E{z[n]z[n  +  k]}  < 
E2{z[n]}  and  pk  <  0  is  achievable. 
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2.1  Case  I  -  Shape  parameter  to  be  an  integer 

The  first  case  considers  an  integer  shape  parameter.  With  this  restriction,  the  method  of  [Armstrong 
&  Griffiths  1991]  first  introduced  the  transform  between  a  correlated  Gaussian  and  correlated 
Gamma.  However,  their  method  only  considered  the  case  of  geometrical  correlation.  In  this 
section,  the  linear  method  is  extended  to  include  any  realisable  correlation  function  but  limited  to 
pk  >  0.  The  criteria  for  realisable  correlations  are  given  in  Subsection  |4.1| 

Since  sequences  with  different  scale  parameters  can  be  obtained  by  re-scaling^J  9  =  1  in  the 
following  discussions.  For  the  special  case  ofz/=l,Z^T(l,l)  and  the  sequence  of  z[n\  is 
found  by  z[n\  =  |x[n]|2,  where  x[n\  is  a  sequence  of  the  Gaussian  process  and  X  ~  CN( 0, 1). 
Therefore, 

(z[n]z[n  +  k])  —  (x[n\x*[n]x[n  +  k\x*[n  +  k])  (3) 

where  the  superscript  *  denotes  complex  conjugate.  Since  x[n\  is  a  complex  Gaussian,  the  ex¬ 
pectation  of  the  right  side  of  ([3])  can  be  found  by  employing  Isserlis’s  theorem  [Michalowicz 
et  al.  2009],  giving, 

(x[n\x*[n\x[n  +  k]x*[n  +  k])  =  (\x\2^  +  \  (x[n]x*[n  +  k}}\2  .  (4) 

Because  ^ \x\ 2^}  —  z2  and  \(x[n\x*[n  +  kj)  |2>  0,  this  results  in  (z[n]z[n  +  k])  >  z2  if  the  real¬ 
isation  of  z[n\  is  through  z[n\  =  |x[n]|2.  Therefore,  regardless  of  the  correlation  of  the  Gaussian 
(x[n\x*[n  +  k]),  pk  >  0  is  the  lowest  bound  of  correlation  coefficient  achievable  for  Z  r(i,i) 
if  its  realisation  is  formed  by  taking  the  intensity  of  correlated  Gaussian  X  ~  CN( 0, 1).  It  can 
be  seen  that  by  taking  the  intensity  of  a  Gaussian  variable  to  form  a  Gamma  variable  inherently 
limits  the  low  bound  of  the  correlation  coefficient  p.  The  correlation  of  Gaussian  is, 

(x[n]x*[n  +  k])  =  d=  ((z[n\z[n  +  k])  —  z 2)1^2  =  =b  (p&var(z))1//2  ,  k  —  0, 1, -  (5) 

Because  p^var^)  >  0,  its  square  root  is  a  non-negative  real  number.  If  x  =  [*[l],®[2],...]r 
is  a  column  vector,  where  the  superscript  T  denotes  transpose,  then 

£{xxH}  =  Mx,  (6) 

where  is  the  covariance  matrix  of  x  and  the  superscript  H  denotes  the  Hermitian  transpose. 
According  to  the  assumptions  of  stationarity  and  symmetry,  M.x  has  a  Toeplitz  structure,  with  el- 
ements  given  by  (|5]),  i.e.,  Mx(n,  k )  =  (p\n-k\)  x  (since  var (z)  =  1  for  Z  ~  T(l,  1)).  According 
to  SIRP,  let, 


x  =  Mi/2u  (7) 

where  U  ~  CA^(0, 1)  is  an  uncorrelated  Gaussian  process.  After  obtaining  the  correlated  x  and 
applying  the  square-law  to  each  element,  the  sequence  (or  the  column  vector)  of 
z  =  [\x  [l]l2,  |a:[2]|2, . .  ,]T  then  has  the  desired  correlation  of, 

E{(z-z)(z-z)t}  =  var(z)Mz  (8) 

3i.e.,  if  Z  ~  r(u,  1),  then  0Z  -  r {v,  0). 


4 


UNCLASSIFIED 


UNCLASSIFIED 


DSTO-TR-2688 


where  M.z  =  Toeplitz (p),  i.e.,  Mz(n,  fc)  =  p\n-k\,  which  is  the  same  as  the  desired  correlation  in 
©• 

The  above  correlated  Gaussian  sequences  can  also  be  obtained  through  Fourier  synthesis 
[Ward,  Tough  &  Watts  2006,  Section  5.4].  The  application  of  a  filter  h  to  a  white  Gaussian  noise 
process  u(t)  gives, 


x(t)  =  I  h(t')u(t')dt' .  (9) 

J  —  oo 

In  the  frequency  domain,  the  above  equation  is  written  as, 

X(f)  =  H(f)U(f )  (10) 

where  £/(/)  is  still  a  Gaussian  process.  On  the  other  hand,  according  to  the  Wiener- Khinchin 
theorem  [Haykin  2007,  Chapter  2],  the  power  spectrum  density  (PSD)  of  the  correlated  signal 

x{t)  is, 


RG(t')e-j27rft' dt’ 


(11) 


where  RG(t)  is  the  auto-correlation  function  of  x(t).  Selecting  H(f)  =  y/Sxx(f)>  the  Fourier 
inversion  of  X(f)  =  )  will  then  result  in  a  correlated  Gaussian  process  x(t)  with  the 

desired  correlation  Rc(t). 

The  advantage  of  the  Fourier  synthesis  against  SIRP  is  that  the  Fourier  transform  in  the  numer¬ 
ical  simulation  can  be  realised  through  the  fast  Fourier  transform  (FFT)  and  hence  significantly 
improves  the  computational  efficiency.  In  this  report,  the  hybrid  method  will  use  the  Fourier  syn¬ 
thesis  instead  of  SIRP  to  generate  correlated  Gaussian  sequences. 

If  v  —  m  >  1,  where  m  is  an  integer,  Lemma  1  can  be  utilised  to  repeatedly  generate  x&, 
k  =  1, . . . ,  m.  The  result  z[n\  =  J2k=i  lx/cM|2  will  then  have  a  r(m,  1)  distribution  with  the 
desired  correlation. 


2.2  Case  II  -  Shape  parameter  to  be  an  integer  plus  0.5 

If  the  shape  parameter  has  a  value  of  v  —  m  +  0.5,  where  mn  >  0  is  an  integer,  the  sequence 
is  repeatedly  generated  for  k  =  1, . . . ,  m  +  1.  The  result  z[n\  =  Ylk=i  lx^N|2  +  Re2  (xm+i  [n]) 
will  then  have  a  T(m  +  0.5, 1)  distribution  with  the  desired  correlation  (according  to  Lemma  1 
and  Properties  1  and  2). 


2.3  Case  III  -  Shape  parameter  to  be  an  arbitrary  number 

If  the  shape  parameter  takes  an  arbitrary  value,  v  =  z/q  +  Az/,  where  =  m  +  0.5/3,  m  >  0  is 
an  integer,  /3  =  0  or  1  and  0  <  Az/  <  0.5.  Two  sequences,  y[n\  and  A y[n\  can  then  be  simulated 
with  respective  shape  values  z/q  and  Az/  so  their  summation  produces  a  sequence  with  the  desired 
shape  and  correlation. 
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The  goal  of  the  hybrid  method  is  to  use  linear  transforms  whenever  possible  to  generate  the 
sequences  of  random  variables.  The  criteria  for  whether  the  linear  method  will  work  is  determined 
by  o  <  rjk  <  1  where  rjk  and  pk  are  related  by, 


r  i  k  =  o 

1  Pk,v/v o  k  ^  0 


(12) 


with  the  correlation  coefficient  rjk  of  the  sequence  y[n\  given  by, 

(y[n]y[n  +  k})  -  y2 


rjk  = 


var  (y) 


k  =  0,1,.... 


(13) 


If  this  relationship  is  valid,  the  hybrid  method  uses  the  Fourier  synthesis  method  to  simulate  the 
desired  sequence.  On  the  other  hand,  it  is  possible  that  this  relationship  will  not  be  satisfied  and  a 
MNLT  technique  must  be  employed.  Each  case  can  be  summarised  by  three  steps: 


Linear  Case 


Step  1:  If  the  relationship  in  (fl2|)  is  satisfied,  generate  a  correlated  Gamma  sequence  of  Y  ~ 
r(i/0,l)  which  has  a  correlation  coefficient  of  rjk- 

Step  2:  Generate  an  uncorrelated  Gamma  sequence  of  AY  ~  r( Az/,  1).  Since  the  sequence 
is  uncorrelated,  the  generation  of  such  a  sequence  is  easy  to  implement  by  simply  using 
either  the  existing  function  in  Matlab  ,  gamrnd,  or  other  MNLT  methods  [Ahrens  & 
Dieter  1974]. 


Step  3:  Combine  the  correlated  sequence  with  r(z/o>  1)  distribution  and  the  uncorrelated  sequence 
with  r(Az/,  1)  distribution  to  produce  z[n\  =  y[n\  +  A y[n\  with  distribution  T(z/,  1). 


To  see  how  the  aggregate  correlation  of  ((y[n\  +  A y[n})  ( y[n  +  k]  +  A y[n  +  k]))  is  equal 
to  the  correlation  of  (z[n\z[n  +  k]},  notice  that  the  sequence  y[n\  is  correlated  whereas 
the  sequence  of  A y[n\  is  uncorrelated  and  the  two  sequences  are  mutually  independent. 
Equation  ([13])  then  becomes, 


9k 


pfevar(z)  +  z2  -  2 A yy  -  y2  -  (AyjnjAyjn  +  k}) 

var  (y) 


(14) 


where  var(z)  uO2,  z  —  vQ,  Ay  —  A i/Q,  y  —  uq9,  var (y)  i/q 92  and 


(Ay[n\Ay[n  +  k]) 


A  v92  +  A  v292  k  =  0 
Av292  k  ^  0 


Inserting  these  values  in  (fl4|)  results  in  ([12]). 


(15) 


MNLT  Case 

Step  1:  If  the  relationship  in  ([12])  is  not  satisfied,  generate  a  correlated  Gamma  sequence  of  Y  ~ 
r(z/0,l)  which  has  a  correlation  coefficient  of  pk- 

Step  2:  Generate  a  correlated  Gamma  sequence  of  AY  ~  T(Az/,  1)  using  the  MNLT  method 
described  in  Section  [2~4l 

Step  3:  Combine  the  correlated  sequence  with  r(z/o?  1)  distribution  and  the  correlated  sequence 
with  r(Az/,  1)  distribution  to  produce  z[n\  =  y[n\  +  A y[n\  with  distribution  T( z/,  1). 
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2.4  A  MNLT  technique  for  shape  parameter  0  <  Ac  <  0.5 

As  mentioned  previously,  not  every  case  of  correlated  Gamma  can  be  tackled  by  the  Linear  method 
and  the  MNLT  must  be  used  instead.  This  subsection  introduces  a  MNLT  technique  specially 
designed  for  tackling  the  case  of0<Az/<0.5. 

The  MNLT  is  a  one-to-one  non-linear  mapping  from  a  known  random  process  (not  necessarily 
Gaussian)  to  another  random  process  with  a  desired  density  function.  Supposing  a  random  process 
y[n\  has  a  known  PDF  py  (y)-  The  MNLT  uses  this  PDF  to  generate  a  random  process  A y[n\  with 
a  desired  PDF  p/yy  (A?/).  This  may  be  realised  by  the  mapping  of, 

coo  coo 

/  pAY(Ay')dAy'  =  /  PY(y')dy'  n  =  1,2,....  (16) 

JAy[n]  Jy[n] 

In  principle,  the  MNLT  method  can  be  used  to  generate  any  non-Gaussian  process.  In  the 
following  discussion  however,  the  method  is  limited  by  using  the  correlated  r(l,  1)  distribution  to 
generate  a  T(Az/,  1)  distribution  with  the  same  correlation. 

Using  the  linear  transform  given  in  the  previous  section,  a  sequence  y[n\  is  first  generated, 
(note  Y  ~  T(l,  1))  with  a  desired  correlation  coefficient  of  pk-  The  one-to-one  non-linear  mapping 
of  (p~6])  is  then  used  to  produce  the  corresponding  A y[n\,  with  a  PDF  of  T(Az/,  1).  This  is  done  by 
equating  their  complementary  distribution  functions  (CDF)  described  in  ([T6|)  giving, 

/  a  \7J  (Au,  A  y[n])  =  exp  (~y[n\) ,  n  =  1,2,...  (17) 

l\Av) 

where  7/ (a,  b )  is  the  upper  incomplete  Gamma  function,  defined  as, 

coo 

7 /(a,  6)=  /  £a_1e_t<iL  (18) 

Jb 

For  a  Gamma  distribution  with  an  integer  shape  parameter,  or  an  integer  plus  0.5,  the  corre¬ 
lation  given  by  ([5])  is  independent  of  the  shape  parameter.  In  another  words,  the  linear  transform 
does  not  change  the  correlation.  According  to  our  simulation  however,  it  was  found  that  the  MNLT 
of  fl7])  does  slightly  alter  the  correlation.  As  a  consequence,  while  the  correlation  of  y[n\  has  the 
desired  values,  after  the  non-linear  transform,  the  correlation  of  A y[n]  will  have  slightly  differ¬ 
ent  values.  A  rigid  mathematical  trace  of  how  the  correlation  is  being  changed  through  the  above 
MNLT  is  difficult  and  hence  an  empirical  formula  is  introduced  to  maintain  the  desired  correlation. 

Consider  the  example  of  generating  a  random  process  A y[n\  using  the  MNLT  of  with 
the  mapping  seeds  of  y[n\.  The  former  has  a  r(Az/,  1)  distribution  with  a  correlation  coefficient 
of  pk ,  and  the  latter  has  a  T(l,  1)  distribution  with  a  correlation  coefficient  of  fa.  The  empirical 
relationship  between  p &  and  fa  is  given  by, 

fa  =  p^Ai/+°'5^  ,  k  =  0, 1, . . .  for  0  <  Av  <  0.5  (19) 

and  is  a  function  of  both  p &  and  Av.  The  transformation  from  p &  to  ^  is  minor  as  shown  in  Figure 
[I]  When  Av  approaches  0.5,  and  p  are  identical,  since  the  transform  from  T(l,  1)  to  T(0.5, 1) 
is  linear. 
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Figure  1:  The  relationship  between  p  and 


2.5  Overview  of  the  Hybrid  Method 

In  summary,  the  proposed  hybrid  method  for  simulating  correlated  Gamma  sequences  primarily 
depends  on  the  shape  parameter  and  may  depend  on  the  correlation  coefficient.  It  is  summarised 
by  the  following  three  cases: 


Case  I  If  v  —  integer,  generate  the  desired  correlated  Gamma  using  the  linear  Fourier  method; 


Case  II  If  v  —  integer  +  0.5,  generate  the  desired  correlated  Gamma  using  the  linear  method; 


Case  III  If  v  is  an  arbitrary  number,  let  n  —  n$  +  An,  where  0  <  Az/  <  0.5  and  z/q  is  an  integer 
or  an  integer  plus  0.5. 


•  If  0  <  Pk^/vo  <  1  is  satisfied,  generate  a  correlated  Gamma  of  shape  z/q  with  the 
correlation  of  (p~3])  following  cases  I  or  II.  Generate  a  second  uncorrelated  Gamma 
sequence  with  shape  An  and  sum  the  two  sequences. 


If  the  check  fails,  generate  a  correlated  Gamma  of  shape  n$  with  the  correlation  of 
pk  following  cases  I  or  II.  Generate  a  second  correlated  sequence  using  the  MNLT 
technique  described  in  Subsection|2.4|to  generate  a  correlated  Gamma  with  shape  An. 
Sum  the  two  sequences  together. 


Therefore,  it  can  be  seen  that  the  generation  of  correlated  Gamma  only  requires  linear  transform  of 
Gaussian  process  for  the  majority  of  cases  to  efficiently  produce  the  correlated  random  sequences. 
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3  Simulation  with  Memoryless  Non-Linear  Transform 

The  concept  of  MNLT  described  in  ([T6|  is  a  non-linear  mapping  from  a  known  random  stochastic 
process  (input)  to  an  unknown  one  (output)  by  equating  their  CDF’s.  If  S  denotes  the  desired 
random  process,  then  its  relation  to  the  Gaussian  process  X  (its  intensity  is  an  Exponential  distri¬ 
bution),  is  given  by, 


where  is  the  PDF  of  S.  The  reason  to  often  choose  a  Gaussian  process  as  the  input,  is  that 
its  correlation  can  be  readily  realised  through  SIRP  or  Fourier  synthesis.  The  difficulty  of  MNLT 
for  correlated  processes  is  that  the  mapping  of  the  correlation  is  also  non-linear.  In  particular,  the 
realisation  of  a  correlated  Gamma  needs  to  first  transfer  the  correlation  of  the  Gamma  into  the 
correlation  of  a  Gaussian.  The  technique  for  manipulating  the  transform  of  the  output  correlation 
back  to  the  input  can  result  in  different  approaches  for  the  MNLT.  In  this  section,  two  approaches 
are  presented.  The  first  is  based  on  using  numerical  integration  and  the  other  employs  inverted 
polynomials  to  map  the  correlation  of  the  output  to  the  correlation  of  the  input. 


3.1  MNLT  Approach  I  -  Numerical  Integration 


This  approach  was  first  described  in  [Tough  &  Ward  1999]  and  expanded  in  more  detail  in  Chap¬ 
ter  5  of  [Ward,  Tough  &  Watts  2006].  The  first  step  of  the  approach  is  to  use  Fourier  synthesis  to 
produce  correlated  sequences  of  random  numbers  with  a  Gaussian  distribution.  The  real  and  imag¬ 
inary  parts  of  the  inverse  Fourier  transform  provide  separate  realisations  of  the  complex  Gaussian 
process.  The  MNLT  which  maps  the  input  Gaussian  random  variables  to  the  output  non-Gaussian 


random  variables  is  via  (20)  and  can  be  written  as, 


£0)  =  Qdist  (erfc  (j^)  /2)  > 


(21) 


where  the  function  Qdist  is  the  complementary  quantile  function  for  the  output  distribution  and 
erfc(-)  is  the  complementary  error  function.  Rapid  evaluation  of  this  function  is  essential  for  prac¬ 
tical  implementation  of  the  method.  For  the  Gamma  distribution,  Qdist  is  the  inverse  incomplete 
Gamma  function.  An  efficient  algorithm  for  this  function  was  developed  and  implemented  in  For¬ 
tran  by  DiDonato  &  Morris  [1986].  The  Fortran  code  and  a  C  translation,  is  available  on  the  web 
as  part  of  the  Cumulative  Distribution  Function  Library  (CDFLIB)  [Venier  &  Serachitopol  2003]. 
The  C  code  has  been  compiled  as  a  Me  x  file  for  Matlab  by  Davidson  [2011]  as  part  of  his  Radar 
Toolbox.  Two  special  cases  are  implemented  with  Matlab  intrinsic  functions:  the  Exponential 
distribution  (y  =  1),  where  Qdist(^)  =  — log(ix),  and  the  Chi-square  distribution  (y  =  1/2), 
where  Qdist(^)  =  erfcinv2^)  and  erfcinv(-)  is  the  inverse  complementary  error  function. 


The  desired  correlation  function  for  the  output  non-Gaussian  process  must  be  mapped  to  the 
correlation  function  for  the  input  Gaussian  process.  The  construction  of  this  mapping  is  described 
in  Section  5.6  of  [Ward,  Tough  &  Watts  2006].  It  requires  the  evaluation  of  integrals  in  the  form 
of  [Tough  &  Ward  1999], 


(((»)((*))  =  /“exp  (-x'2)  ffpMQfc,  <fa',  (22) 
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where  Ro(t)  =  (x(0)x(t))  /  (x2)  is  the  correlation  coefficient  of  x(t)  (note  x  is  real  and  X  ~ 
N( 0, 1))  and  Hp(-)  is  the  Hermite  polynomial  of  order  p.  The  Hermite  polynomial  coefficients 
are  obtained  with  the  function  HermitePoly  written  by  David  Terr,  found  on  Matlab  Central 
[Terr  2004].  The  integrals  are  evaluated  with  12  point  Gauss-Hermite  integration,  using  the  func¬ 
tion  Gauss  Hermite,  which  is  a  MATLAB  translation  of  the  Fortran  function  hermite  .  f  90 
by  Miller  [2010]. 

The  advantage  of  this  MNLT  is  that  Gamma  sequences  with  p(t)  <  0  can  be  simulated,  though 
the  correlation  is  still  be  bounded  to  the  condition  given  in  Subsection  |4.1|  According  to  Chapter 
5  of  [Ward,  Tough  &  Watts  2006],  for  the  case  of  v  =  1,  ([22])  reduces  to, 

<£(0)£(*)>  ~  1  +  0.816 RG(t)  +  0.177i?G(t)2  +  0.0067JRG(t)3 

+0.00013.RG(t)4  +  0.000017.RG(t)5.  1  ; 

This  means  that  a  Gamma  process  with  a  correlation  p{t)  <  0  can  be  simulated  by  choosing  a 
proper  Rc(t)  <  0. 


3.2  MNLT  Approach  II  -  Polynomial  Approximation 

When  the  MNLT  converts  a  correlated  Gaussian  process  into  a  new  process  with  desired  marginal 
distributions,  the  new  process  is  also  correlated,  but  not  through  the  same  auto-covariance  function 
as  the  Gaussian  process.  It  is  in  fact  related  to  the  Gaussians  auto-covariance  by  a  non-linear  map¬ 
ping.  In  a  practical  situation,  a  process  with  a  given  auto-covariance  function  would  be  specified. 
It  is  shown  that  by  using  an  appropriate  inversion  method,  the  correlation  of  the  Gaussian  can  be 
obtained  with  the  desired  correlation  for  the  new  process  [Weinberg  &  Gunn  2011a]. 

As  mentioned,  the  MNLT  method  is  an  extension  of  the  CDF  inversion  method  with  which 
a  single  random  variable  can  be  generated.  For  a  given  continuous  random  variable,  inverting  its 
CDF,  and  evaluating  it  on  uniformly  distributed  random  numbers  between  0  and  1,  will  generate 
sequences  distributed  from  the  original  random  variable.  This  can  be  used  to  very  easily  produce 
independent  samples.  To  extend  this  to  the  generation  of  correlated  random  samples,  it  is  required 
that  if  X\  and  X2  are  two  random  variables,  with  distribution  functions  Fx1  and  Fx2  respectively, 
then  the  random  variable 

f£(FXi(Xi))  ~X2.  (24) 

The  validity  of  this  result  is  discussed  in  [Weinberg  &  Gunn  2011a].  Equation  ([24])  is  the  key 
to  generating  correlated  realisations  of  a  random  process.  This  can  be  seen  by  defining  a  function 
=  Fxl  (Fx1  (x))  and  observing  that  if  x  is  a  realisation  of  random  variable  X\,  then  £(x)  is 
a  realisation  of  random  variable  X 2. 

Consider  a  wide  sense  stationary  stochastic  process  ((t)  evolving  over  time  t ,  with  correlation 
function  R^(r)  =  E(£(0 )C(r)),  and  each  £(£)  ~  X\.  This  can  be  transformed  to  produce  a  new 
stochastic  process  6{t)  —  £(£(£)),  which  will  be  distributed  according  to  X2  pointwise.  However, 
this  new  process  will  have  a  correlation  function  Rq(t)  =  k(R^{t)),  where  k(-)  is  a  non-linear 
function. 

The  advantage  of  this  is  that  any  random  variable  X2  can  be  generated  from  any  X\.  If  X\ 
is  used  to  generate  a  correlated  Gaussian  sequence  which  is  then  transformed  into  a  sequence 
of  random  variables  with  the  desired  marginal  statistical  distribution  using  the  function  £.  The 
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resulting  sequence  will  then  consist  of  dependent  realisations.  Due  to  the  non-linear  nature  of  £, 
the  output  sequence  will  not  necessarily  have  the  same  correlation  properties  as  the  Gaussian  input 
process.  Consequently,  the  relationship  between  input  and  output  auto-covariance  functions  is 
examined  in  detail  in  [Tough  &  Ward  1999].  In  particular,  if  the  Gaussian  process  auto-covariance 
is  denoted  Rq,  and  the  desired  correlated  process  has  auto-covariance  R0 ut,  then  ([22])  can  be 
written  as, 


Rout(t)  — 


oo 

E 

p= 0 


(flc(*))p, 

2 Pp\ 


"pi 


with  the  coefficients  ap  given  by 


1  poo 

=  —=  /  e~x  Hp(x)£(V2x)dx. 

V^r  J _ oq 


(25) 


(26) 


The  relationship  ([25])  would  be  more  useful  if  it  could  be  inverted,  to  specify  the  Gaussian  auto¬ 
covariance  necessary  to  generate  a  desired  auto-covariance  in  the  output  process.  Some  simple 
examples  where  this  relationship  can  be  inverted  are  considered  in  [Tough  &  Ward  1999]. 


This  issue  is  addressed  in  detail  in  the  report  [Weinberg  &  Gunn  2011a].  In  particular,  a 
polynomial  approximation  to  ([25])  is  derived,  which  is  inverted  and  scaled,  and  consequently  used 
to  determine  an  appropriate  Gaussian  input  process.  As  an  example,  four  terms  were  manually 
calculated,  resulting  in  a  polynomial  approximation  of, 


Rg(p)  ~  57  (-Rout(p)  -  bo) 

-ftCRout(p)-M2 

-  (l  -  f )  (^out (p)  -  bo)3  (27) 

-(|-^  +  f)  (^out (p)-bo)\ 

where  bp  =  ap/(2pp\). 

This  has  been  found  to  work  well  in  a  number  of  cases  considered.  Furthermore,  the  validity 
of  this  polynomial  approximation  has  been  examined  in  [Weinberg  &  Gunn  2011  b].  In  terms 
of  practical  simulation  of  desired  correlated  processes,  it  has  been  reported  that  the  success  of 
this  approach  can  vary  greatly.  According  to  our  simulation,  it  has  been  found  that  the  approach 
produces  poor  results  for  the  correlated  Gamma  with  small  shape  parameters,  v  <  0.3.  This  may 
be  due  to  the  size  of  the  polynomial  used  for  the  approximation.  Further  investigation  is  required 
for  this  case. 
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4  Examples  and  Comparisons 

A  number  of  different  methods  have  been  described  for  simulating  correlated  Gamma  random 
variables.  Correlation,  however,  has  its  own  physical  and  mathematical  constraints  which  have  not 
been  discussed.  Section |47T| presents  two  constraints  to  ensure  that  a  given  correlation  function  is 
mathematically  (and  possibly  physically)  realisable.  Two  examples  of  applying  the  hybrid  method 
are  then  presented  in  Section  |4.2|  A  performance  comparison  of  the  methods  is  then  given  in 
Section|43j  Finally,  a  description  of  how  the  hybrid  method  can  be  used  to  simulate  K-distributed 
data  is  given  in  Section  [4~4| 


4.1  Realisable  Correlation  Functions 

This  section  presents  two  constraints  which  are  required  to  ensure  that  a  correlation  functions  is 
realisable. 

Constraint  1:  As  shown  by  ([2])  the  correlation  coefficient  of  a  correlated  Gamma  sequence  z[n\ 
satisfies  po  =  1  and  \pk\  <  1  for  k  >  0. 

Constraint  2:  The  correlation  matrix  given  by  ([8])  is  positive  semi-definite.  The  proof  is 
given  below. 

Under  the  stationary  and  symmetric  assumptions,  the  covariance  matrix  of  Gamma  dis¬ 
tributed  random  vector  z  =  [z[  1], . . . ,  z[iV]]T  is  given  by, 

X  =  E  j(z  —  z)  (z  —  z)H^  =  var(z)M^.  (28) 

For  an  arbitrary  constant  (non-random)  vector  q  =  [qi, . . . ,  gA]T,  the  quadratic  form  is, 

qIi'Eq  =  E  {(qH(z- z))  (( z-z)Hq )}  =  E  {\qH  (z  -  z)f}  >0.  (29) 

Therefore  the  covariance  matrix  X,  or  equivalently  the  correlation  matrix  (because  of 
var (z)  >  0),  is  positive  semi-definite. 

These  two  constraints  are  the  characteristics  of  Gamma  distributed  sequences  (wide-sense 
stationary  sequences)  and  independent  of  simulation  methods.  If  the  hybrid  method  is  used  to 
generate  correlated  gamma  sequences,  the  aforementioned  additional  condition  of  pk  >  0  is  also 
required. 

The  validity  of  the  correlation  using  the  above  conditions  must  therefore  be  checked  when 
simulating  a  correlated  Gamma  process  with  a  given  set  of  correlation  coefficients.  Any  correlation 
not  satisfying  the  above  constraints  cannot  be  realised  accurately  by  the  above  methods.  For 
example,  consider  the  correlation  coefficient, 

Pk  —  |cos(0.l7rfc)|  |fc|  =  0, 1, -  (30) 

A  simple  check  indicates  that  is  not  positive  semi-definite  and  hence  the  desired  correlation 
cannot  be  precisely  generated  by  any  of  the  above  simulation  methods.  However,  for  a  correlation 
coefficient  of, 

Pk  =  cos(0.1257rfc)e-fc//1°,  |fc|  =  0, 1, . . . ,  (31) 
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the  corresponding  M.z  is  positive  definite  and  hence  the  desired  correlated  correlation  can  be 
simulated  by  the  MNLT  methods.  However,  due  to  the  condition  of  pk  >  0  not  being  satisfied,  the 
proposed  hybrid  method  is  unable  to  produce  the  desired  correlation. 


4.2  Examples  with  the  Hybrid  Method 


To  demonstrate  the  hybrid  method,  two  numerical  examples  are  now  presented.  The  first  has  a 
correlation  coefficient  of, 

pk  =  [0.7  +  0.3cos(0.127rfc)]  e~k/12,  k  =  0,1,2,....  (32) 

The  desired  shape  parameter  is  v  =  3.7  which  gives  v o  =  3.5  and  Av  =  0.2.  A  simple  check 
shows  the  condition  of  ijk  —  pk  v/vq  <  1  is  satisfied.  Therefore,  the  linear  transform  was  used  to 
generate  a  correlated  dataset  that  has  a  T(3.5, 1)  distribution  with  a  correlation  coefficient  of  pk- 
The  second  uncorrelated  dataset  having  a  T(0.2, 1)  distribution  was  simply  produced  by  calling 
the  Matlab  function  gamrnd.  The  combination  of  the  two  will  then  give  a  T( 3.7, 1)  distribution 
with  the  desired  correlation  coefficient  of  pk.  Even  though  typical  radar  systems  may  only  collect 
a  few  hundred  range  bins  at  a  time,  a  sequence  containing  106  samples  was  generated  to  ensure 
there  are  sufficient  samples  for  confirming  the  actual  correlation.  The  comparison  between  the 
desired  and  simulated  correlations,  probability  and  cumulative  densities  are  shown  in  Figure  [2j  It 
can  be  seen  that  the  desired  correlation  of  the  Gamma  has  been  precisely  achieved.  The  second 


Simulated  values  (dB) 


Figure  2:  Comparison  between  the  desired  and  simulated  correlation  coefficients  for  Example  1: 
( — )  exact,  (—)  simulated:  (top)  one-sided  correlation  coefficients ;  ( middle )  PDF’s  of  the  Gamma 
and  ( bottom )  CDF’s  of  the  Gamma. 
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example  has  a  correlation  coefficient  of, 

pk  =  e-k/w,  k  =  0,1,2, ... . 


(33) 


In  this  case,  the  desired  shape  parameter  is  v  —  0.7  which  gives  vq  =  0.5  and  Av  =  0.2.  A 
simple  check  shows  the  condition  of  %  =  pk^/u o  <  1  is  not  satisfied.  Therefore,  the  linear 
transform  is  used  to  generate  a  correlated  dataset  that  has  a  Y (0.5, 1)  distribution  and  the  MNLT 


method  described  in  Subsection  2.4  is  used  to  generate  a  correlated  dataset  that  has  a  Y (0.2, 1) 
distribution.  Both  of  them  have  a  correlation  coefficient  of  p^.  The  combination  gives  a  dataset 
that  has  a  T( 0.7, 1)  distribution  and  the  desired  correlation.  The  comparison  between  the  desired 
and  simulated  correlation  coefficients,  probability  and  cumulative  densities  are  shown  in  Figure[3] 
The  results  again  show  that  the  desired  correlation  of  the  Gamma  has  been  precisely  achieved. 


Figure  3:  Comparison  between  the  desired  and  simulated  correlation  coefficients  for  Example  2: 
( — ),  exact  (-■-)  simulated:  (top)  one-sided  correlation  coefficients ;  ( middle )  PDF’s  of  the  Gamma 
and  ( bottom )  CDF’s  of  the  Gamma. 

4.3  Comparison  of  Algorithms 

In  this  section,  comparisons  among  the  hybrid  method,  numerically  integrated  MNLT  approach  as 
well  as  the  polynomial  MNLT  implementation  are  performed  in  terms  of  accuracy  and  timeliness. 
In  particular,  the  shape  parameter  is  varied  from  0.1  to  10  in  steps  of  0.1  and  the  mean  is  set  to  1. 
104  samples  are  generated  for  each  realisation  and  the  result  is  averaged  100  times.  It  was  noted 
that  the  polynomial  MNLT  method  does  not  work  efficiently  below  a  shape  of  0.3. 
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The  first  comparison  looked  at  the  accuracy  of  the  measured  auto-correlation  function  by  com¬ 
paring  the  simulated  with  the  desired.  It  was  found  that  each  of  the  three  methods  are  very  accurate 
with  the  maximum  RMS  error  over  the  100  iterations  being  less  than  0.002. 

The  second  comparison  looks  at  the  relative  error  in  the  shape  parameter  and  scale  by  sub¬ 
tracting  the  maximum  likelihood  estimate  from  the  true  value  and  dividing  that  result  by  the  true 
value.  Figure  [4]  shows  the  results  with  the  shape  in  the  top  plot  and  scale  in  the  bottom.  The  shape 
parameter  estimate  is  quite  good  for  each  method  with  a  small  relative  error.  However,  there  is  a 
slight  increase  in  the  relative  error  for  the  polynomial  MNLT  method  when  the  shape  is  below  0.7. 
The  scale  estimate  is  also  very  accurate  for  all  methods  when  the  shape  is  above  1 .  For  low  shape 
parameters,  however,  the  hybrid  method  is  more  accurate. 

The  third  comparison  looks  at  the  second,  third  and  fourth  order  moments  of  the  simulated 
data.  These  are  also  known  as  the  variance,  skewness  and  kurtosis.  Figure  [5]  shows  that  all  three 
methods  produce  accurate  moments  with  the  numerically  integrated  MNLT  method  being  less 
accurate  at  low  shape  parameters. 

An  exact  expression  for  the  algorithm  complexity  is  difficult,  as  the  exact  implementation 
details  are  not  known  for  all  algorithms.  Instead,  the  final  comparison  looks  at  the  mean  run  times 
for  the  three  methods.  Figure  [6]  shows  that  the  polynomial  MNLT  method  is  the  slowest  by  an 
order  of  magnitude,  while  the  performance  of  the  hybrid  method  varies  depending  on  the  shape 
value.  The  hybrid  method  is  the  fastest  for  nearly  all  shape  values  and  approaches  the  numerically 
integrated  MNLT  method  as  the  shape  increases.  This  improvement  is  significant  considering  that 
the  inverse  incomplete  Gamma  function  in  the  latter  method  was  coded  in  C  in  order  to  improve 
computational  efficiency,  while  the  hybrid  method  only  uses  builtin  Matlab  functions. 

In  summary,  the  hybrid  method  is  slightly  more  accurate  at  low  shape  values  than  the  other 
MNLT  methods  and  has  a  superior  mean  run  time.  The  polynomial  MNLT  method  is  less  accurate 
with  low  shape  parameters  and  takes  a  significantly  longer  time  to  run. 
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Figure  4:  Parameter  estimate  error:  ( — )  hybrid ,  (-X-)  MNLT  poly.,  (—)  MNLT  numerical. 
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Shape 

Figure  5:  Moment  error:  ( — )  hybrid,  (-X-)  MNLT poly.,  (—)  MNLT numerical. 


Shape 

Figure  6:  Mean  run  time  comparison:  ( — )  hybrid,  (-X-)  MNLT  poly.,  (—)  MNLT  numerical. 
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4.4  K-distribution  Example 

As  a  final  example,  two  range-time  (pulse)  intensity  images  of  correlated  K-distributed  sea-clutter 
are  shown  in  Figure  [7]  The  temporal  correlation,  i.e.  the  correlation  among  the  fast- varying 
Gaussian  component  with  respect  to  pulse,  is  specified  by  ([33]),  while  the  spatial  correlation,  i.e. 
the  correlation  among  the  slowly-varying  Gamma  component  with  respect  to  range  bin  is  specified 
by  <[32]). 
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(b)  v=3.7 


Figure  7:  Range-time  (pulse)  intensity  maps  of  correlated  K-distributed  sea-clutter.  Both  have  the 
same  correlation  but  different  shape  parameters.  Normalised  intensity  values  >0.5  are  shown  as 
white  in  the  images. 
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5  Conclusion 

This  report  has  systematically  studied  how  to  generate  correlated  Gamma  processes  to  support  the 
simulation  of  sea-clutter  for  many  DSTO  projects.  In  particular,  a  hybrid  method  was  proposed  for 
generating  the  desired  correlated  Gamma  samples.  The  method  uses  only  linear  transformations 
for  most  values  of  the  shape  parameter  and  correlation,  with  a  non-linear  transform  required  in 
some  cases.  The  difficulty  of  the  non-linear  transform  is  the  mapping  of  the  correlation  of  the  de¬ 
sired  output  to  the  correlation  of  the  input,  often  a  Gaussian  process.  Two  MNLT  implementations 
were  presented  for  comparison,  one  based  on  numerical  integration  and  the  other  a  polynomial 
approximation.  Performance  of  the  three  methods  has  been  evaluated,  focusing  on  the  accuracy  of 
the  desired  correlated  process  and  the  timeliness  for  the  realisation.  It  was  found  that  the  hybrid 
method  is  the  slightly  more  accurate  for  small  shape  parameters.  It  is  also  the  most  computation¬ 
ally  efficient  and  hence  the  fastest  for  nearly  all  shape  values.  This  improvement  is  significant 
considering  that  the  numerically  integrated  MNLT  method  has  the  inverse  incomplete  Gamma 
function  coded  in  C  in  order  to  improve  computational  efficiency,  while  the  hybrid  method  only 
uses  builtin  Matlab  functions.  The  drawback  of  the  hybrid  method  is,  however,  that  it  can  only 
handle  positive  correlations  whilst  the  two  MNLT  methods  are  capable  of  handing  both  positive 
and  negative  correlations. 
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Appendix  A  Temporal  and  Spatial  Correlations  of 
Compound  K-Distributed  Clutter 

The  compound  K-distribution  assumes  that  the  data  consists  of  a  fast-varying  component  modu¬ 
lated  by  a  slowly-varying  component.  The  fast- varying  component  commonly  refers  to  the  speckle 
component  that  is  a  Gaussian  random  process  with  zero  mean  and  unit  variance.  The  slowly- 
varying  component  also  describes  the  underlying  mean  or  texture  with  intensity  modelled  by  a 
Gamma  distribution.  Care  must  be  taken  when  estimating  the  correlations  of  these  components 
as  the  slowly-varying  component  may  not  always  remain  constant  when  estimating  the  temporal 
correlation. 

A.l  Temporal  Correlation 

The  temporal  correlation  is  present  between  data  samples  collected  by  a  pulse  train  in  a  coherent 
processing  interval  (CPI).  Often  during  the  CPI,  it  can  be  assumed  that  the  slowly-varying  com¬ 
ponent  remains  unchanged,  i.e.,  the  slowly  varying  component  is  fully  correlated  (this  is  how  the 
compound  K-distributed  clutter  got  its  name Q  The  covariance  matrix  of  the  data  is  written  as, 

Wi  —  E  {xx^}  ,  (Al) 

where  expectation  is  with  respect  to  the  pulse  for  the  temporal  correlation,  and 
x  =  [^[0],  •  •  •  ?  x[N  —  1]]T  is  an  N  x  1  vector  collected  by  N  pulses.  Each  measurement  may 
be  further  written  as  a  product  of  fast-varying  and  slowly-varying  components,  according  to  the 
model  assumption, 


x[n\  =  y/rxf[n\,  (A2) 

where  Xf[n\  is  the  fast- varying  component  which  is  a  complex  Gaussian  variable  and  r  is  the 
underlying  mean  which  is  constant  in  a  CPI.  Because  the  two  components  are  independent, 

M  —  E  {xx^}  =  E  {r}  E  {xjx^}  /iMj,  (A3) 

where  M f  =  E  jxjx^j  is  the  covariance  matrix  of  the  fast- varying  component  and  /i  is  the 
clutter  mean.  Therefore,  the  temporal  correlation  can  be  estimated  using  data  samples  and  the 
correlation  of  the  fast-varying  component  is  just  the  covariance  matrix  normalised  by  its  mean. 
Under  the  assumption  of  wide  sense  stationary,  M  f  has  a  Toeplitz  structure  of, 


= 


1 

Pi 

Pn-i 


Pi  •  •  •  Pn-i 

1  •••  : 

:  '  • .  pi 

•••  P\  1 


where  p^  —  \ E  {x[n\x*[n  +  k]}  ,  n,  k  =  0, 1, . . . ,  N  —  1. 


(A4) 


4With  reference  to  a  maritime  search  radar,  there  is  high  superposition  of  antenna  footprints  with  respect  to  succes¬ 
sive  pulses  in  a  CPI,  so  the  texture  of  sea-clutter  can  be  assumed  to  be  completely  correlated. 
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A.2  Spatial  Correlation  (Correlation  in  Range) 

In  order  to  find  the  correlation  of  the  slowly  varying  component,  the  data  needs  to  be  manipu¬ 
lated  in  the  intensity  (or  amplitude)  domain  rather  than  the  complex  (I/Q)  domain.  Consider  the 
correlation  between  range  bin  n  and  n  +  k,  k  ^  (Q 

E  |  y/r\n\x  f  [n]  VTln  +  k]x*f[n  +  k]  j  —  E  j  \Jr[n}r[n  +  k\  j  E  {xf[n]x^[n  +  £;]}  .  (A5) 


If  however,  the  interval  between  bin  n  and  bin  n  +  k  is  greater  than  the  radar  range  resolution, 
E  ^Xf[n]x*f[n  +  k\  j  =  0.  Therefore,  even  if  E  |^/r[n]r[n  +  k]  j  /  0,  its  value  is  not  measur¬ 
able  using  the  I/Q  data.  The  correlation  therefore  has  to  be  found  using  the  intensity  (or  amplitude) 
data.  The  correlation  coefficients  of  texture,  r[n\  and  intensity,  z[n\  =  x2[n]  are  given  by, 


Pk 


E  {r[n]r[n  +  k]}  —  E 2  {r} 
var(r) 


&  =  0,1,..., 


(A6) 


with 


Xk 


E  {z[n\z[n  +  k]}  —  E 2  {z} 
var  (z) 


k  =  0,1,..., 


(A7) 


Xo  =  Po  =  I-  (A8) 

The  next  step  is  to  find  the  relationship  between  \k  and  pk  f°r  k  ^  0.  Denoting  z  =  r\xf\2  =  rzf, 
gives 


E  {z[n]z[n  +  k}}  —  E  {r[n]r[n  +  k]zf[n]zf[n  +  k]}  (A9) 

and  since  the  fast- varying  and  slowly  varying  components  are  independent, 

E  {r[n\r[n  +  k]zf[n]zf[n  +  k]}  =  E  {r[n]r[n  +  k]}  E  {zf[n\zf[n  +  k]}  .  (A10) 


For  k  ^  0,  the  value  of  E  {zf  [n\zf  [n  +  k] }  can  be  calculated  using  Isserlis’  Theorem  [Michalowicz 
et  al.  2009],  giving 

E  {zf[n\zf[n  +  k]}  —  E  |x/[n]x*  [n]xf[n  +  k\x^[n  +  fc]| 

=  1  +  E  $yXf[n]Xf[n  +  fc]|  E  +  /c]j  for  k  ^  0.  (All) 

+E  | Xf[n]x*f[n  +  k]j  E  |x^[n]xj[n  +  k\  | 


If  the  interval  between  bin  n  and  bin  n  +  k  is  greater  than  the  radar  range  resolution,  the  last  two 
items  of  (|A11|)  become  zero  and  E  {zf[n\zf[n  +  k}}  =  1,  simplifying  the  above  correlation  to 


E  {r[n]r[n  +  k]}  =  E  {z[n]z[n  +  k]}  for  k  /  0.  (A12) 


5In  order  to  simplify  symbols,  the  same  index  notation  is  used  to  represent  either  temporal  series  (the  index  refers  to 
pulse  numbers)  or  spatial  series  (the  index  refers  to  range  bin  number).  There  should  be  no  confusion  under  the  context. 
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Combining 


and  (|A12|)  and  noticing  E{r}  =  E{z}  =  fi,  then  gives 
var(r) 


var(z)  ’ 


k  =  1,2. 


(A13) 


Equations  flA8[)  and  (|A13|)  indicate  that  once  the  correlation  of  z  or  r  is  known,  the  correlation  of 
the  other  can  be  determined.  The  texture  r  is  Gamma  distributed,  and  the  intensity  z  is  single-look 
or  multi-look  K-distributed,  with  variances  given  respectively  as, 


var(r)  =  fi2/u, 


var  (z) 


v  +  N  +  1 

ivv 


(A  14) 
(A15) 


where  N  is  the  number  of  multi-looks.  Therefore,  due  to  the  effect  of  the  fast- varying  component 
that  is  uncorrelated  and  randomly  varying,  the  correlation  of  the  intensity  z  is  generally  weaker 
than  the  correlation  of  the  texture  r.  Only  if  the  number  of  multi-looks  reaches  infinity,  will  the 
fluctuations  of  the  fast- varying  component  disappear  (averaged  to  its  mean  value  for  every  range 
bin)  and  the  two  correlations  become  identical. 


The  final  relationship  is  therefore  found  by  inserting  (|A14|)  and  (|A15|)  into  (|A13|),  giving, 


N 

Xk  ~  Pku  +  N  +  V 


k  =  1,2.... 


(A16) 


The  correctness  of  the  above  derivation  is  confirmed  by  Monte  Carlo  simulation  where  the  slowly 
varying  Gamma  component  has  a  shape  parameter  of  v  —  1.2  and  a  correlation  coefficient  of, 

pk  =  [0.7  +  0.3cos(0.127rfc)]e~fc/12,  A:  =  0,1,....  (A17) 


The  multi-look  K  distributed  data  was  then  generated  by  modulating  the  uncorrelated  multi-look 
fast- varying  Gaussian  component  with  the  correlated  Gamma.  The  correlation  coefficient  for  the 
simulated  multi-look  K  data  was  regressed  and  compared  to  the  theoretical  value.  Figure  |A.2 
shows  a  comparison  between  the  ideal  and  theoretical  correlation  for  the  single  look  case  and 
with  16  looks.  It  can  be  seen  that  the  simulated  results  match  the  theory  with  the  multi-look 
averaging  reducing  the  variance  of  fast- varying  component.  As  a  result,  <  \x\ 2  >  1  and  hence 

the  correlation  approaches  the  correlation  of  texture.  For  the  single-look  case,  oscillation  of  the 
uncorrelated  fast-component  greatly  reduces  the  overall  correlation. 

For  over-sampled  range  data,  when  the  range  interval  between  bin  n  and  bin  n  +  k  (k  /  0)  is 
smaller  than  the  range  resolution,  the  fast- varying  component  is  also  correlated  in  range,  resulting 
in  E  {zf[n]zf[n  +  k]}  >  1  (see  flAll[)).  The  correlation  of  the  intensity  z  will  then  be  jointly 
contributed  by  the  fast- varying  component  and  the  slowly-varying  component. 
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(a)  Correlation  of  gamma  (b)  correlation  of  K 


Figure  Al:  The  desired  correlation  of  the  (a)  Gamma  component  (16  looks )  and  (b)  the  correlation 
of  the  multi-look  K  data  for:  (~o~,  -*-)  1  look  and  (-o-,  -*-)  16  looks. 
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